/*
 * =====================================================================================
 *
 *       Filename:  point_detection.h
 *
 *    Description:  point detection
 *
 *        Version:  1.0
 *        Created:  2009年06月08日 18时58分52秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */

#ifndef POINT_DETECTION_H
#define POINT_DETECTION_H
#include "point.h"
#include "process.h"
#include <vector>
typedef std::vector<NCpoint<int> > NCpointVector;

/**
 * |rx2  rxry|
 * |rxry ry2 |
 */
NCpointVector harris_corner (const NCmatrix<double> &mat, const NCmatrix<double> &kernel,double k,double thr)
{
	int i,j,hi,hj,row=mat.row(), column=mat.column(),hrow=kernel.row(),hcolumn=kernel.column();
	double w,rx,ry,rx2,ry2,rxry;
	NCmatrix<double> gradx = gradientx(mat), grady = gradienty(mat);
	NCpointVector pointVec;
	for( i=0; i<row; i++ )
	{
		for( j=0; j<column; j++ )
		{
			rx2=0;ry2=0;rxry=0;
			for( hi=0; hi<hrow; hi++)
			{
				for ( hj=0; hj<hcolumn; hj++ )
				{
					if ( 
							(i-(hrow-1)/2+hi)>0    && 
							(j-(hcolumn-1)/2+hj)>0 && 
							(i-(hrow-1)/2+hi)<row    && 
							(j-hcolumn/2+hj)<column 
					   )
					{
						w = kernel[hi][j];
						rx = gradx[i-(hrow-1)/2+hi][j-(hcolumn-1)/2+hj];
						ry = grady[i-(hrow-1)/2+hi][j-(hcolumn-1)/2+hj];
						rx2 += w*rx*rx;
						ry2 += w*ry*ry;
						rxry += w*rx*ry;
					}
				}
			}
			if( (rx2*ry2-rxry*rxry-k*(rx2+ry2)*(rx2+ry2)) > thr )
			{
				pointVec.push_back(NCpoint<int>(i,j));
			}

		}
	}
	return pointVec;
}
#endif
